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Abstract 

A first core is a first hydrostatic object formed in the course of dynamical con- 
traction of a molecular cloud core. Since the inflow pattern changes drastically before 
and after the first core formation, it is regarded as a milestone in the star formation 
process. In order to identify the first core from a mapping observation, the features 
expected for the first core are studied for CS rotation transitions at radio wavelengths. 
The non-LTE radiation transfer is calculated for the results of radiation magnetohy- 
drodynamical simulations of the contraction of the magnetized molecular cloud core 
in rotation (Tomida et al. 2010a). We use the Monte-Carlo method to solve the non- 
LTE radiation transfer in a nested grid hierarchy. In the first core phase, an outflow 
arises from the vicinity of the first core due to the twisted magnetic field amplified 
by the rotation motion of the contracting gas disk. The disk and outflow system has 
several characteristic observational features: (i) relatively opaque lines indicate asym- 
metry in the emission lines in which the blue side is stronger than the red side (an 
infall signature of the envelope); (ii) in the edge-on view, the disk has a signature of 
simultaneous rotation and infall, i.e., the integrated intensity of the approaching side 
is brighter than that of the receding side and the gradient in the intensity-weighted 
velocity is larger in the approaching side; (iii) the observed outflow indicates rotation 
around the rotation axis. The size of the outflow gives the approximate age after 
the first core is formed, since the outflow is not expected for the earlier runaway 
isothermal collapse phase. 

Key words: stars: formation — ISM: jets and outflows — radiation transfer - 
ISM: lines and bands — methods: numerical 
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1. Introduction 

The star formation process continues to attract much attention. In the process, "the 
first core" is regarded as a key to understanding the star formation process, since it is a first 
hydrostatic object created in the course of forming a new star. The first core was theoretically 
predicted by Larson (1969) as an object supported by gas pressure after the thermal radiation 
from the dust, which is a major coolant in interstellar molecular gas with density n ~ 10 4 — 
10 10 H2cm~ 3 , becomes inefficient at a high density of n ^ 10 10 ~ 11 H2 cm" 3 . The high density gas 
increases in temperature due to gravitational compression and a hydrostatic balance is finally 
achieved between the self-gravity and the pressure gradient. This forms a hydrostatic object, 
the first core. 

The dynamical evolution is completely different before and after the first core formation 
(Whitworth & Summers 1985; Tomisaka 1998). Before the first core formation, the contraction 
is well expressed by the so-called Larson (1969)-Penston (1969) (LP) similarity solution. After 
the core formation, the solution resembles Shu (1977) 's inside-out collapse solution but is well 
expressed by another self-similar solution regarded as an analytic continuation of the LP solu- 
tion, as given in Whitworth & Summers (1985). Magnetohydrodynamical simulations confirm 
the transition at the first core formation even in a model with a magnetic field and rotation 
motion (Tomisaka 1998; Machida, Tomisaka, & Matsumoto 2004). After the core formation, 
the magnetic field lines, which are essentially poloidal before the first core formation, begin to 
be distorted by the rotation motion and the toroidal component is amplified. This transfers its 
angular momentum from the rotating disk to the gas just outside the first core and accelerates 
the gas, which explains the molecular outflows observed accompanying protostars. Therefore, 
the formation of the first core is a milestone in the process of star formation. 

Although more than 40 years have passed since Larson's (1969) prediction, the first core 
has not been observed. Recently however, a possible detection has been reported in a dense 
core IRS2E in L1448 (Chen et al. 2010). An object has been discovered which is not visible with 
Spitzer (3.6 — 70/im) but is a weak mm or sub-mm continuum source. Saigo & Tomisaka (2011) 
calculated the expected spectral energy distribution (SED) of dust thermal emissions for their 
hydrodynamic model of the first core (Saigo, Tomisaka, & Matsumoto 2008). They obtained a 
luminosity of ~ 0.02 L and effective temperature of ~ a few 10 K. These expected features of 
the first core, i.e., low luminosity and cold SED, are consistent with the observation by Chen 
et al. (2010). However, this object seems more evolved than the first core since it is associated 
with a high- velocity outflow (~ 25km s" 1 ). Another candidate is Per-Bolo 58, a dense core in 
Perseus, found in Spitzer 70 fim map (Enoch et al. 2010). Although this has also a similar SED 
to that expected from the first core, detection of 24 /im seems to indicate the object is older 
than the first core and the envelope is being to be cleared. Therefore, further observations are 
required to confirm whether these objects are first cores. Saigo, Tomisaka, & Matsumoto (2008) 
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predicted that a first core forms a flat disk due to rotation and that a dynamical instability 
develops in such a fast rotating disk (Durisen et al. 1986). This brings us an idea that a 
fine structure (such as spriral arms) should be observed in first core disks in high resolution 
observations with such as ALMA (Tomida et al. 2010b; Saigo & Tomisaka 2011). 

Another signature of a first core is an accompanying molecular outflow (Tomisaka 1998; 
2002). The molecular outflow is driven by the magnetic Lorentz force working just outside the 
first core (Tomisaka 2002; Banerjee & Pudritz 2006; Machida, Inutsuka & Matsumoto 2007; 
Commergon et al. 2010; Tomida et al. 2010a). In this paper, we investigate the appearance 
of the molecular outflow in its early phase, which should allow identification of a first core. 
This is done by post-processing the numerical results of radiation MHD simulations. We call 
this procedure "observational visualization." It is a method of visualizing numerical results in 
order to understand the undergoing physics but also emphasizes the observational expectation 
from the simulation. We have already shown the expected polarization observation of the 
thermal emission from dust grains in the molecular outflow (Tomisaka 2011). Comparison 
with observation, this enables us to conclude whether the molecular outflow is directly driven 
magnetically or is indirectly made by the entrainment mechanism from the interstellar jets. 

The plan of the paper is as follows: In section 2 we describe the model and numerical 
method. A non-LTE line transfer calculation is made for previously obtained data from ra- 
diative MHD simulations of the contraction of a rotating magnetized cloud. The model of the 
initial molecular cloud is described. The method for the non-LTE radiative transfer calculations 
for the interstellar molecular lines is also described in this section. In section 3, we show the re- 
sults of the observational visualization. Section 4 is devoted to a discussion of the asymmetrical 
distribution with respect to the vertical axis and a comparison with the isothermal model. 

2. Model and Method 

Tomida et al. (2010a) (hereafter Paper I) have calculated the evolution of a rotating 
magnetized cloud using radiation magnetohydrodynamical (RMHD) simulations. In Paper I, 
we assumed a cloud in a hydrostatic balance with a central density of p c = 1.0 x 10~ 19 gcm~ 3 or 
no = 1.6 x 10 4 H 2 cm -3 , a uniform rotation rate of u = 0.1/tfj = 1.5 x 10 _14 rad s _1 , and relatively 
weak uniform magnetic field of B z = 1.1 /iG. The angular momentum vector and the magnetic 
field direction are assumed to both be in the ^-direction. 

In Paper I, we employed the nested grid technique, in which the global structure is cov- 
ered with a coarse grid while the central spatially fine structure with high density is calculated 
with a fine grid. We numbered the grid level from L = (the coarsest) to L = 17 (the finest). 
The structures of the respective grid levels are self-similar. All the levels of the grids are co- 
centered, i.e., the centers of all the grids are placed at the same point (see Fig. 2). The number 
of grid cells of each level was chosen to be 64 3 and the size of cell A^ of the L-th level is given 
by that of the L = level A as A L = 2~ L A . 
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Before the second collapse phase n < 10 15 H2 cm -3 , the main coolant of the molecular gas 
is the dust thermal emission, rather than the gaseous line cooling emissions. RMHD simulations 
of Paper I included the radiation transfer of the thermal emissions using the gray flux-limited 
diffusion (FLD) approximation. We apply here the radiation transfer calculation of gaseous 
line emissions to the snapshot data of density, kinetic temperature, and velocity distributions 
obtained in Paper I. Thus, we employ a so-called post-process to calculate the molecular line 
emissions. Figure 1 shows the density (a)-(c) and temperature (d) distributions obtained in 
Paper I. The ages of these two stages are t = 3.852 x 10 5 yr (a and b: just before the first core 
formation) and t = 3.858 x 10 5 yr (c and d: r = 6.45 x 10 2 yr after the first core formation). That 
is, Figures 1(a) and (b) represent the prestellar isothermal collapse phase, while Figures 1(c) 
and (d) correspond to the protostellar first core phase. At the stage of Figure 1(c) and (d), the 
first core has grown to the mass of ~ 0.04M Q . From the first core, an outflow with 1 — 2 kms" 1 
is ejected, of which the velocity is similar to the Kepler speed of the launching point (Kudoh 
& Shibata 1997). 

2.1. N on- LTE Radiation Transfer 

We calculate the level population of the rotation transitions of a number of abundant 
molecules in molecular clouds and solve the non-LTE radiative transfer problem for the rotation 
transitions of the molecules. Denoting the number density of the J level (energy level E(J)) 
as nj, we can write the balance equation as 

nj ^2 Rjj, = ^2 nj,Rj,j ( J = 0, 1, • ■ ■ , J max ), (1) 
j'^j j'^j 

where Rjj' represents the transition probability from J to J 1 as 

f = Ajj, + Bjj,J v jj, + nCjj, for J > J', 
JJ ' \ = Bjj,J v jj, + nCjj, for J < J', 

where Ajji and Bjj, represent Einstein's coefficients, the former being the coefficient for spon- 
taneous emission and the latter the coefficient for absorption (J < J') and induced emis- 
sion (J > J'). Cjji is the collisional transition rate from J to J' for collisions with H2 
molecules whose density is denoted by n. The average intensity of radiation with a frequency of 
v = [E(J') — E(J)]/h is written as J v jj>, where h is the Plank constant. We take into account 
the energy levels from J = to J = J max . Although J max should be taken as large as possible 
for completeness, in our simulations of cloud collapse and outflow we take J max = 10. 

To obtain the average intensity J U) we have to solve the radiation transfer equation for 
the specific intensity I u : 

-^- = -K v I v + e y , (3) 

where k u and e„ represent the monochromatic volume absorption coefficient and monochromatic 
volume emissivity. Total absorption coefficient and emissivity for the transition between J and 



J' are related to the level populations rij as 



(tot) _ hv 
JJ ' ~ 4nAu 



{nj,Bj,j-njBj.j,) (J>f), (4) 



^ = ^ njAjJ ' {J>JI) - (5) 
Here the Doppler width of the line due to the thermal and turbulent motion of molecules is 

expressed as 

Au = (6) 

c 

using the microturbulence parameter a. Assuming a line profile function <f) of 



v - VQJJ' 



where vqjji = [E(J) — E(J')]/h, the absorption coefficient and emissivity for v are written using 
the total absorption coefficient and emissivity for specific transition J — > J' as 

KuJj' = (j)j.ji{v)n { jf) (8) 

£yjj> =<pjj>{v)t { j°r ■ (9) 

Thus, solving equations (1) and (3) is a nonlocal problem, or in other words, physical states of 
the molecules in different places are coupled with each other by the radiation. 

We solve this non-LTE radiative transfer problem by a Monte Carlo method similar to 
van der Tak (2000); Hogerheijde & van der Tak (2000). The formalism is as follows: The region 
is divided into uniform Cartesian cells of A^ id . In randomly chosen directions, rays are ejected 
from the center of each cell. Along the ray, the transfer equation (3) is solved from the outer 
boundary to the cell center. Before we finish integrating equation (3) for all rays from each 
cell, the level population, n"j ld \ and thus the emissivity and the absorption coefficients are 
both fixed. The number of rays per cell is chosen to be = 100 in this calculation. After 
completing the integration, we can obtain the average intensity J for each cell by 

J JVRay 

Jujj> = f Iujj'dQ = —— J2 lyJJ'- ( 10 ) 

J 7V Ray N=1 

This gives new level populations, nj iew ' ) , consistent with the average intensity J for use in 
equation (1). This cycle is repeated until the assumed level populations n^° ld ^ and the new 
populations ?7,j Lew ' ) converge. We assume the convergence criterion that a relative deviation 
of the populations is smaller than e n = 10 -8 , \rij — n j iew ' 1 1 jn j° ld ^ < e n . The procedure is 
tested by a comparison with a test problem done by Juvela (1979) which has been used for the 
visualization of interstellar turbulence driven by supernovae (Wada & Tomisaka 2005; Yamada 
et al. 2007). 

Molecular data, such as the Einstein's A and B coefficients and the C coefficients in 
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equation (2), are taken from the Leiden Atomic and Molecular Database (LAMDA) 1 (Schoier 
et al. 2005). The abundance of CS relative to H2 is chosen to be Xcs = 4 x 10~ 9 . Frequencies 
in laboratory frame of the respective transitions are shown in Table 1. The one- dimensional 
random velocity composed of the thermal plus turbulent contributions is chosen to be a = 
100 ms" 1 of equation (6) 2 . 

2.2. Radiative Transfer on Nested Grid 

In this subsection we present the method to solve the non-LTE radiative transfer on the 
nested grid hierarchy. When the cloud is in a vacuum, the outer boundary condition should 
be that the inwardly directed intensity is equal to that of the Cosmic Microwave Background 
(CMB; Tb = 2.7K). In reality, the interstellar molecular cloud is immersed in the interstellar 
radiation field. However, since we do not have sufficient information about the interstellar 
radiation field, we adopt the vacuum boundary condition in this paper. In order to reduce the 
artificial effect of this boundary condition, we place the boundary as far as we can from the 
inner part which we are interested in by using the nested grid technique. As shown in §2, the 
nested grid technique uses a special grid hierarchy as the grid spacing increases with departing 
from the center. This enables us to place the outer boundary far from the center compared 
with the uniform grid spacing. In the nested grid hierarchy, we assume the outermost grid 
L = is in a vacuum filled with the CMB. The outer boundary of the inner nested grids (L > 1) 
should be taken from the intensity obtained from the coarser outer grids (see Fig. 2). In this 
case, we calculate the level populations of the molecules in each nested grid step-by-step from 
the coarser level to the inner finer level as shown in Appendix. 

Also, in the 3D non-LTE simulations, the number of cells in each level of the nested 
grid is chosen to be equal to that of the RMHD simulation, Nf D = 64 3 . The grid sizes of 
these simulations are chosen to be the same, A^. The cell size is taken as A^ = 1150AU/2 L . 
Thus, in the model shown in Figure 1 (L = 9) the size is equal to A9 = 2.25 AU. We use the 
Cartesian coordinate system (x,y,z) in calculating non-LTE radiative transfer, which is also 
used in RMHD simulation (Fig. 3 left) . 

In this paper, we consider axisymmetric objects. To reproduce the observations, we 
choose a special geometry in which (1) the direction of the line of sight (LOS) is on the x-z 
plane and (2) in the figures the horizontal direction coincides with the y-axis (see Fig. 3). Under 
this specification, the LOS is expressed as n = (sin^O^os^), where the direction of the LOS is 
specified only with the angle measured from the z-axis, 8, since the structure is axisymmetric. 
If we express the directions of the observational grid as e.\ (horizontal) and e.2 (vertical) (that 

1 http://www.strw.leidenuniv.nl/~moldata/ 

2 We changed a from 100ms -1 to 300ms" 1 . Although obtained line-width and line-shape depend on a, results 
are qualitatively the same, because the global velocity gradient has more significant effects on the optical 
depth than the local line width. 
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is, any observational point is expressed with two integers i and j as x = Ax,(iei + je.2)), the 
unit vectors are given as e.\ — (0,1,0) and — [e g — (e z ■ n)n]/\e z — (e z • n)n| (see Fig. 3). The 
case for 9 = 0° is a pole-on view of the disk and that for 9 = 90° is edge-on. 

2.3. Efficiency of Nested Grid 

To clarify the efficiency of the nested grid we compare results obtained with and without 
the nested grid radiation transfer for the model of Figure 1. The integrated intensity distribution 
(color) 



where V and Tb represent the velocity along the LOS and the corresponding brightness tem- 
perature, are shown for CS J = 1-0 (left panels: a and c) and J = 8-7 (right panels: b and 
d) lines in Figure 4. The upper panels (a and b) are obtained without the nested grid method 
for L = 9 density, temperature, and velocity distributions, to which we apply the CMB bound- 
ary condition on the outer boundary of L = 9. In contrast, the lower panels (c and d) show 
the result with the nested grid technique. That is, the outer boundary is properly placed 
on the L = grid and the intensity at the outer boundary of L = 9 is obtained consistently. 
Comparing (a) and (c), we can see that the integrated intensity of the J = 1-0 line seems to 
be very centrally peaked in the calculation without the nested grid method. Both the gradient 
in the integrated intensity and the velocity gradient V(V) are overestimated in the calcula- 
tion without the nested grid method. This shows that not accounting for the envelope outside 
L = 9 has a strong impact on the intensity distribution for a relatively optical thick case of the 
J = 1-0 line. On the other hand, the difference in the integrated intensities between (b) and 
(d) is relatively small compared with that between (a) and (c). This is due to the low optical 
thickness in the envelope for this high excitation line J = 8-7, since such emissions occur only 
in the central part of the disk and the J = 7 level is not excited in the outer low-temperature 
envelope. This seems to explain the similarity in (6) and (d). The comparison indicates that 
it is important to use proper boundary conditions in the non-LTE calculations. That is, even 
if the excitation temperature T ex is obtained accurately (for example, LTE is established) for 
a grid level L without the nested grid technique, the optical thickness r outside the level L 
cannot be ignored in estimating the line intensity. In other words, the intensity of the optically 
thick lines is qualitatively inaccurate if we ignore the outer grids. That is to say, even in such 
a case, without applying the nested grid method, a proper results for the non-LTE radiation 
transfer is not obtained. 

Since we assume 9 = 60° here, the obtained structures seen in (a) and (b) must be affected 
by the fact that the LOSs passing the central part and those passing uppermost and lowermost 




(11) 



and the intensity weighted velocity (contours) 




(12) 
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parts have different path lengths in the L = 9 level. However, we place the outer boundary far 
from the target level of the grid, L = 9, by the nested grid technique. Thus, this effect of the 
artificial boundary is removed in (c) and (d). 

3. Results 

3.1. Spectral Change between Prestellar and Protostellar Phases 

In Figure 5, we show the CS position-to-position spectra of J = 2-1 and J = 7- 6 for 
the prestellar isothermal collapse phase [(a) and (b)] and for the protostellar first core phase 
[(c) and (d)}. The area covered by the spectra is the same as the grid of L = 9 shown in Figure 
1. These spectra are for the edge-on view of the disk (9 = 90°). The difference between the 
stages before and after the core formation is clearly seen in this 140- AU scale spectra. In the 
isothermal collapse phase, since the gas is essentially isothermal with a temperature of Xk — 10 
K, the brightness temperature is equal to Tb — 7-8K. In the first core phase, the peak brightness 
temperature reaches Tb ~ 30K for J = 2-1 and Tb ~ 25K for J = 7-6. In this scale, the first 
core with a high temperature contributes to the spectra. (In Figure 1 (d), gas with Tk ^ 20K 
extends to r ^ 30 AU from the center.) 

The line width also indicates a clear difference. The line width of the first core phase (full- 
width half- maximum ~ 2kms _1 ) is apparently larger than that of the isothermal collapse phase 
(~ lkms -1 ). Figures 5 (a) and (b) indicate a two-peak signature. The blue- and red-shifted 
components represent, respectively, the approaching and departing sides of the contracting gas. 
Since, in the contracting gas, the near side gas with a positive recession velocity experiences 
self-absorption due to the foreground contracting envelope, the signature in which the blue 
peak (far side) is stronger than the red peak (near side) indicates that the gas is contracting 
(Zhou et al. 1993). This shows that the brightness of the CS J = 2-1 line includes contributions 
from the infalling gas not only in the isothermal collapse phase but also in the first core phase. 
The difference between ( a) and ( c) shows that the infall speed is accelerated from the prestellar 
to the protostellar phase. Outflow, which begins to appear in the first core phase, is seen as 
a "low-velocity" component in the positions (±l,+4), (±1,-4) in panels (c) and (d). This 
component can be traced near the center between (+1,-2) and (+l,+2). 

3.2. Effect of the Viewing Angle 

In Figure 6, we show the integrated intensity distributions and intensity-weighted mean 
velocity (V) for the CS J = 2-1 and J = 7-6 lines for the prestellar isothermal collapse phase 
of Figure 1 (a). In this subsection we consider the effect of the viewing angle 9 of LOS. 

From the models of CS J = 2-1, we can see that the spatial variation of the integrated 
intensity in CS J = 2-1 is as small as ~ 20%. This is due to the fact that the line is relatively 
optically thick and the velocity gradient is not so large in this phase. Therefore, the integrated 
intensity does not follow the density distribution (as was already shown from the similarity in 
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the CS J = 2-1 spectra in Fig.5). 

Panel (a) shows the pole-on view with 6 = 0° and panel (k) shows the edge-on view with 
9 = 90°. The direction of the minor axis of the integrated intensity coincides with the rotation 
axis. The CS J = 7-6 line has a low integrated intensity from 6.6 to 7.8Kkms _1 , compared with 
J = 2-1 (18-19 Kkms -1 ). However, an antisymmetric velocity pattern, which features global 
rotation motion, is well traced in CS J = 7-6 but not in J = 2—1. The velocity pattern specific 
to the rotation motion is erased in the CS J = 2-1 emission owing to its optical thickness. 

Figure 7 is the same as Figure 6 but for the protostellar first core phase. The model 
with 9 = 0° (pole-on: panels (a) and (b)) exhibits a ring-like enhancement of intensity. This 
enhancement corresponds to the outflow lobe seen in Figure 1 (c) and (d). It should be noted 
that overall (V) is negative when the system is viewed pole-on. However, since the ring has a 
more negative recession velocity than the rest of the field, this negative velocity originates from 
the outflowing gas. 

Viewing from 9 = 30° [(c) and (d)}, the feature becomes more complicated. Namely, in 
addition to the ring-like feature a bar-like structure appears near the center in the integrated 
intensity distribution. In the negative velocity region, an island of positive velocity region 
appears at — 25AU < x < 25AU and — 30AU < V ~ — 5AU. There are two possible explanations 
for this positive velocity region: (l)infalling gas entering the disk with positive recession velocity 
absorbs the emission coming from the hotter interior. If such gas is evacuated by the outflow, 
the blue-shift emission is strengthened. However, if this is the case, a similar signature must 
appear in the pole-on model (a and b). Thus we consider the second explanation, (2) a positive 
recession velocity from the accelerated outflow in the far side. This positive velocity region 
moves to the right with 9 and is merged with the signature of the rotating disk (the left-hand 
side has a negative velocity and the right-hand side has a positive velocity) for 9 > 60° [(g)-(j)]- 
Finally, the edge-on view [(h) and (/)] indicates a fat disk-like appearance in the integrated 
intensity distribution which is approximately symmetric against the rotation axis. It should be 
noted that the gradient in (V) is larger than that of the isothermal collapse phase 3 . 

However, this distribution has a slight asymmetry in the brightness. Globally, the (V) 
distribution is antisymmetric against the rotation axis, which apparently indicates that the 
disk is rotating. However, this also exhibits a deviation from the antisymmetric distribution. A 
larger portion of the gas seems to have negative velocity. This asymmetry against the rotation 
axis is due to the fact that the gas has both rotation and infall motions, which is seen more 
clearly in §§ 3.3 and 4.2. 

3.3. Rotation and Infall Motion 

To investigate the origin of the asymmetric distributions, we show the integrated inten- 
sity and intensity- weighted average velocity for the CS J = 2-1 (a) and J = 7-6 (6) lines in 

3 The step of the isovelocity contour was chosen to be 0.025 kms -1 in Fig. 6 but is 0.05 kins -1 in this figure. 
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Figure 8. To compare these, we made the same plots for a rotation model (c and d) in which we 
preserved the rotation motion and artificially removed the infall velocity and for an infall model 
( e and f) in which we preserved the infall motion and artificially removed the rotation velocity. 
The rotation and inflow velocities are calculated using the tangential unit vector = e r x e z 
as 



respectively 

As expected, the rotation and infall models have completely different integrated intensity 
and the intensity- weighted velocity distributions. In the infall model (c and d), the integrated 
intensity distribution exhibits a compact core in a fat spheroid. The observed first moment 
velocity is symmetric against the z-axis. In contrast, in the rotation model (e and f), a flat disk 
is seen in the integrated intensity distribution and an outflow extends to z < 70AU in the first 
moment velocity, which indicates an antisymmetric pattern. Thus, an accurate model in which 
gas rotates and simultaneously inflows has no symmetry except against the x-axis (y = 0). This 
asymmetry can be qualitatively understood from the fact that the superposition of symmetric 
and antisymmetric velocities has no symmetry. 

We should also note from Figure 8 (a and b) that the gas to the left of the plot is relatively 
brighter than the gas on the right (a and b). However, the infall and rotation models have 
symmetric integrated intensity distributions against the rotation axis. This can be understood 
as being due to the effect of rotating inflow gas. We discuss the reason in more detail in section 
4.2. 

4. Discussion 

4-1. Comparison with Isothermal Model 

In Figure 9, we compare the radiation model [(a) and (b)] with the isothermal model 
[(c) and (d)], in which the temperature distribution is artificially replaced with the isothermal 
distribution of kinetic temperature Tk = 10 K. Panels (a) and (c) show the CS J = 2-1 line 
and (b) and (d) show the CS J = 7-6 line. The figure indicates that the isothermal model 
has extremely low brightness even in the protostellar first core phase (the CS J = 2-1 line 
has a peak of 25.5 Kkms" 1 while the radiative model predicts 89.6 Kkms -1 . The CS J = 7-6 
line has a peak of ~ 10 Kkms" 1 in the isothermal model while the radiative model predicts 
~ ISOKkms" 1 ). Since the peak of the integrated intensity is at ~ ^Kkms" 1 for J = 2-1 and 
~6Kkms~ 1 for J = 7-6 in the isothermal collapse phase, this isothermal model predicts a similar 
brightness as for the isothermal collapse phase. This indicates that if we ignore the radiation 
transfer in determining the kinetic temperature the emission is significantly underestimated in 
the protostellar phase, at least at this scale, which is not unexpected. 





(13) 
(14) 
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Another difference between the two models is that the velocity gradient in the disk is 
significantly underestimated in the isothermal model. That is, the contour of (V) ~ O^kms" 1 
runs at |x| ~ 55 AU in (a) but \x\ ~ 40 AU in (c). This is a natural outcome of low temperature 
for the isothermal model Tk = 10 K. At this temperature the molecules are essentially in 
the ground state. The effect of the optical thickness does not play an important role in the 
isothermal model, in which the integrated intensity follows the column density distribution and 
the isovelocity seems to follow the mass average velocity. This clearly shows that both the gas 
distribution and the gas dynamics are likely to be misinterpreted if we ignore the radiation 
transfer. 

4-2. Asymmetry against the Rotation Axis 

In this subsection we consider why the left side of the disk and outflow is observed to 
be brighter than the right side. We assume a gas with both contraction and rotation motions, 
in which the infall v± n and rotation speeds depend only on the radius r as v ; n (r) and v^r). 
Consider two LOSs A-B and A'-B', symmetric with respect to the center O, as shown in 
Figure 10. Points A and A' are tangential points of the LOSs. Asymmetry arises from the 
configuration that emission from the gas with higher T ex near the tangential point is absorbed 
by the foreground cooler gas. The velocity difference between the emitter and the absorber has 
an essential role in the self-absorption. The LOS velocities at points A and A' are respectively 
v^Tq) = v^q and — ity(ro) = — v^q, where r represents the distance between A (and A') and O. 
Those at points B and B' are equal to v</,(r) cosa + t>i n (r) sina and — v^(r) cosa + v in (r) sina, 
respectively. Here L AOB is denoted by a and r represents the distance between B (and B') 
and O. The relative velocity of point B observed from point A becomes 

Av AB = v^(r) cosa - + v- m (r) sina, (15) 

while that of A'B' is equal to 

Ava>b> = -i^(r)cosa + f^o + 'Win(^)sina. (16) 

Since v in (r) since > 0, Avab > Ava'b 1 and |Auab| > lAi^^l, when ity(r) cos a — v^q > 0. When 
v<j,(r) cosa — v^o < 0, Avab < Ava'B' and |Avab| < \ Ava>b'\- This is valid irrespective of a. For 
rigid body rotation, the rotation speed satisfies v$(r) cosct — v^o — 0. Since the rotation law of 
this kind of object is between the Kepler rotation oc r -1 / 2 and the rigid body rotation oc r +1 , 
the inequality v^{r) cos a — < is satisfied for r > r . In this case, the magnitude of the 
relative velocity for A'B', |Aua'b'|, is larger than that for AB, |At>AB|,that is, | Ava'b> \ > \ Avab\- 
This means that the approaching side (A'B') has a larger velocity gradient than the receding 
side (AB), when inflow and rotation motions coexist. 

Since the excitation temperature decreases with radius, the emissions from the inner 
radii (A and A') are more or less absorbed by the intervening foreground gas (B and B'). 
Absorption depends not only on the excitation temperature but also the velocity difference 
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between the emitter and absorber if we consider a line with a relatively large optical depth that 
shows a "self-absorption" feature. A larger velocity gradient leads to a less efficient absorption. 
In this case, since |Aua'S'| > |Auab|, self- absorption due to the foreground gas is less efficient 
for the LOS A'B'. As a result, when the gas has both infalling motion and rotation motion, the 
approaching side (A'B') is brighter than the receding side (AB). This explains the asymmetry 
around the rotation axis. That is, the rotation and simultaneous contraction motions of a gas 
leads to an approaching-side-enhanced asymmetry, if the excitation temperature decreases with 
radius. 

4-3. Observability of the First Core Stage 

To find the first core, we have to know how it can be observed. In Figures 6 and 7, we 
have shown how the observational features change after the first core is formed. In the face-on 
view, we have found that a bright spot and a ring are typical for the first core phase. The former 
indicates a first core and the latter indicates molecular outflow driven magnetically. Since the 
second core is believed to have a much brighter intensity owing to its deeper gravitational 
potential, objects found with the intensity calculated here ~ 50Kkms _1 are regarded as being 
candidates for the first core. In the edge-on view, the first core phase has a typical appearance 
of a disk and magnetically driven molecular outflow, both of which have strong blue-intense 
asymmetry. This indicates strong evidence that the rotation and infall motions coexist. In 
Figure 9, the lowest level of the integrated intensity represents the outflow at the 70- AU scale. 
The height of the molecular outflow is approximately proportional to the age after the first 
core is formed, since the molecular outflow began to be accelerated just after the first core 
formation. If we measure the length of the molecular outflow as < 100 AU, the central core 
must be a first core. 

5. Summary 

We have calculated the non-LTE radiation transfer for CS rotation transitions, which 
gives an expectation of the observation features of the first core. A Monte-Carlo code is devel- 
oped to solve the non-LTE radiation transfer based on the nested grid hierarchy. Incorporating 
the nested grid hierarchy enables us to calculate the case in which a large optical depth is 
expected for the envelope in the foreground. Viewing from the rotation axis, the line width 
allows identification of the transition between before and after the first core formation. The 
spectra of the disk show blue-intense asymmetry, which indicates the inflow. The molecular 
outflow launched from the vicinity of the first core appears as a 50-AU-scale ring in the pole-on 
view and is observed as a rotating funnel in the edge-on view. The rotating and inflowing disk 
has a characteristic feature which is seen in neither a purely rotating disk nor a purely inflowing 
disk: The approaching side is brighter than the receding side in the integrated intensity and 
the velocity gradient is larger in the approaching side than in the receding side. This enables 
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us to identify the rotating inflowing disk. This asymmetry is also seen in the molecular outflow 
gas. 

This work was supported in part by JSPS Grant-in-Aid for Scientific Research (A) 
21244021 in the fiscal year 2009-2010. K. Tomida was supported by the Research Fellowship 
from JSPS for Young Scientists. Numerical computations were carried out in part on NEC SX-9 
and Cray XT4 at the Center for Computational Astrophysics, CfCA, of National Astronomical 
Observatory of Japan. 
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Appendix. Non-LTE Radiation Transfer Calculation on Nested Grid Hierarchy 

In this Appendix, we describe the numerical method to solve the non-LTE radiation 
transfer calculation on the nested grid hierarchy. The nested grid uses a number of grids whose 
spatial resolutions are different. The finer grid covers a central region and the coarser one 
covers the whole cloud (Fig. 2 a). Each level of grid contains the same number of grid points. 
The boundary condition for the coarsest (level L = 0) grid is given as the cloud is immersed 
in a vacuum which is filled with the CMB. On the contrary, the boundary condition for the 
grids contained in the L = grid (L > 0) is determined by the strength of intensity I u at the 
boundary which is calculated using the emissivity and absorption coefficient of the outer grids. 

1. Solve L = grid with the CMB condition at the outer boundary. Non-LTE radiative 
transfer calculation gives distributions of emissivity and absorption coefficient obtained in 
L = level (Fig.26). 

2. Solve L = 1 grid. 

(a) Generate rays passing L = 1 grid. 

(b) Integrate radiative transfer equation along the rays for the part covered by L = grid 
but outside L = 1 grid. Record the intensity at the outer boundary of L = 1. 

(c) By iteration, obtain the equilibrium solution for L = 1 grid with the boundary con- 
dition obtained in the previous step (2-b), giving the distributions of emissivity and 
absorption coefficient in L = 1 level (Fig.2c). 

3. Solve L = 2 grid. 

(a) Generate rays passing L = 2 grid. 

(b) Integrate radiative transfer equation along the rays inside L = grid but outside L = 2 
grid using emissivity and absorption coefficients obtained in both L = and L = 1 
grids. When a point is covered both by L = and L = 1 level grids, the emissivity 
and absorption coefficients for the finer level (L = 1) are taken according to the basic 
idea of the nested grid. Record the intensity at the outer boundary of L = 2. 

(c) With iteration, obtain the equilibrium solution for L = 2 grid with the boundary 
condition obtained in the previous step (3-b). This gives the distributions of emissivity 
and absorption coefficient in L = 2 level (Fig. 2d). 

4. When the target grid level is reached, we stop this procedure. 
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Transition CS Frequency (GHz) ALMA Band 

J =1-0 48.9909549 

J = 2- 1 97.9809533 Band 3 

J = 7 -6 342.8828503 Band 7 

J = 8-7 391.8468898 Band 8 

Table 1. Frequencies in laboratory frame of the transitions appeared in this paper. 
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Fig. 1. Density structure of prestellar isothermal collapse phase at 2200-AU (a: Level 5) 
and 140- AU (6: Level 9) scales. The contour levels are H2 number density n = 10 6,5 cm~ 3 , 

10 cm , 10 cm , The central density reaches n c ~ 10 10 H 2 cm at this time. In 

this phase, gas is essentially isothermal with ~ 10K, which is not shown in this figure. In 
(c) we show the density distribution after the first core formation. The kinetic tempera- 
ture is plotted in (d) with contour levels T = 12. 5K, 15K, 17. 5K, The size of pan- 
els (c) and (d) is 140 AU x 140 AU (Level 9). The velocity field is shown by arrows. 
We plot a velocity vector of (v x ,v z ) = (lkms -1 ,0) in the lower-left corner for comparison. 
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Fig. 2. Method of solving non-LTE radiative transfer problem on the nested grid, (a) Geometry 
of the nested grid, where the L = 1 grid covers the central 1/4 of the L = grid, the L = 2 grid 
covers the central 1/4 of the L = 1 grid, and so on. Although 2-dimensional cross-cut is shown 
in this figure, actual radiative transfer problem is solved in 3-dimensional geometry, (b) In this 
method, we first solve the non-LTE problem of L = with a CMB outer boundary condition. 
This gives the distributions of the absorption coefficient and emissivity in the L = grid, which 
enables us to calculate the intensity at the outer boundary of L = 1. (c) We solve the non-LTE 
problem of L = 1 with the intensity at the boundary obtained in (b) . This gives the distributions 
of absorption coefficient and emissivity in L = 1 grid, which enables us to calculate the intensity 
at the outer boundary of L = 2. (d) We solve the non-LTE problem of L = 2 with the intensity at 
the boundary obtained in (c). This procedure continues until we have reached the target level. 
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Non-LTE and RMHD 
simulation grid 



Observation 
grid 



Fig. 3. Relationship between the grid used in the non-LTE and RMHD simulations (left: nested 
grid) and the observation grid (right). Observations were made by integrating along the normal vec- 
tor n. That is, the direction of the observation is specified by n = (sin#, 0, cos^). The unit vectors 
which specifies the observation grid are given as e\ = (0,1,0) and &2 = [e z — (e z -n)n]/|e z — (e z ■ n)n\. 
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Fig. 4. Comparison of expected integrated intensities with and without applying the nested grid method. 
The integrated intensities based on the non-LTE radiation transfer for the model of Fig.l are shown by 
the false color. These are the CS J = 1-0 ([a] and [c]) and J = 8-7 ([b] and [d]) transitions viewed 
from 6 — 60°. The contour lines represent the intensity- weighted mean velocity (V). The upper pan- 
els ([a] and [b]) indicate the case where the extracted L = 9 level data are integrated without apply- 
ing the nested grid method. In the lower panels ([c] and [d]) the radiation transfer is calculated using 
all the data for levels < L < 9. Color bars near the upper- left corner represent the levels of inte- 
grated intensities in Kkms -1 . The x- and y-axes represent e\ and e-i axes in Fig. 3, respectively. 
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Fig. 5. Spectra observed from edge-on 6* = 90° in the L = 9 level. Each spectrum is delineated by the ob- 
served position. The upper panels (a and b) represent the model before core formation (isothermal phase) 
and the lower panels (c and d) represent the model after the core formation (first core phase). The horizon- 
tal axis represents the recession velocity (— 1.5kms _1 < V < l.Skms" 1 ). The left panels (a and c) are for CS 
J = 2-1 and the right ones (6 and d) are for CS J = 7-6. The maximum of the vertical axis is taken to be 40 
K. The positions of the spectra are labeled using the x and y distance from the center. That is, the spec- 
trum in the upper- right corner is (4,4) and the upper-left one in the central four spectra is labelled (—1,1). 
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(a) (b) (c) (d) 




(e) (f) (g) (h) 




(i) 0) (k) (1) 




Fig. 6. CS J = 2-1 and J = 7-6 integrated intensity (false color) and intensity-weighted mean velocity 
(V) for different LOSs (contours) based on the data of Fig. 1(a). The panels (a), (c), (e), (g), (i), and (fc) 
are results for the CS J = 2-1 lines and panels (b), (d), (f), (h), (j), and (t) are results for the CS J = 7 6 
lines. The viewing angles 9 for the respective panels are (a) and (b): 6 = 0° (pole-on), (c) and (d): 9 = 30°, 
(e) and (f): (9 = 45°, (#) and (h): 9 = 60°, (i) and (j): 9 = 80°, and (fc) and (I): 9 = 90° (edge-on). The levels 
of the integrated intensity are shown in the color bar in the upper- left corner and the unit is Kkms -1 . 
The solid and dashed contour lines of (V) represent positive (Okms -1 < (V) < 0.5 kms -1 ) and negative 
(— O.Skms" 1 < (V) < Okms" 1 ) velocities, respectively. The step of the contour is chosen to be 0.025kms _1 . 
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(e) 
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Fig. 7. The same as Fig. 6 but for the protostellar first core phase (Fig.l (c) and (d)). The velocity 
range is Okms^ 1 < (V") < lkms -1 for positive (solid-line contour) velocity and — lkms -1 < (V) < Okms -1 
for negative (dashed-line contour) velocity. The step of the contour is chosen to be 0.05 kms -1 . 



23 



(b) 



data_bif_T/Lev09 

77.3700 
J=2-1, th=90 



data_b1f_T/Lev09 

66.8900 
J = 7-6, th = 





(c) 



(d) 



60 



20 



data_b1f_inf/Lev09 

91.9100 
J— 2-1, th=90 



data_b1f_inf/Lev09 





















/ " 

y 




V 




-1 




/ 

/ 






X 

\ 






/ 

/ 

/ 


/ 

It 
\ 


ms-g 

~c c 
- c V 




2\' f X N 




20 


\ 














\ 


\ 






< . ° / 
/ 




-20 




K \ 

N 

"N 






"~ ~ / 

/ 

/ 














/ 




-60 


-60 


-40 


-20 





20 40 60 







c . '-'if- = XX\ 4 
< ' c X^ - X c <^ - 



-60 -40 -20 



20 40 60 



x(AU) 



x(AU) 



(f) 



data_b1f_rot/Lev09 
66.0800 



data_b1f_rot/LevG9 
73.6500 





Fig. 8. Comparison between models in which only the inflow motion is assumed [(c) and (d)) and that 
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